Superfluid to Mott-insulator transition of hardcore bosons in a superlattice 



Itay HcrEl and Marcos RigofI 
Department of Physics, Georgetown University, Washington, DC 20057, USA 
(Dated: October 20, 2009) 

We study the superfluid to Mott-insulator (SF-MI) transition of hardcore bosons in commensurate 
superlattices in two and three dimensions. We focus on the special case where the superlattice has 
period two and the system is at half filling. We obtain numerical results by using the stochastic 
series expansion (SSE) algorithm, and compute various properties of the system, such as the ground- 
state energy, the density of bosons in the zero-momentum mode, the superfluid density, and the 
compressibility. We employ finite-size scaling to extrapolate the thermodynamic limit, and find the 
critical points of the phase transition. We also explore the extent to which several approximate 
solutions such as mean-field theory, with and without spin-wave corrections, can help one gain 
analytical insight into the behavior of the system in the vicinity of the phase transition. 

PACS numbers: 64.70.Tg, 03.75.Lm, 02.70.Ss, 67.85.-d 
Keywords: superfluidity, Mott-insulator, hardcore bosons 



I. INTRODUCTION 

Recent developments in the field of ultracold Bose 
gases have opened a new promising avenue of theoret- 
ical and experimental research in the study of the phases 
of quantum matter. A gas of bosonic atoms in an opti- 
cal trap has been reversibly tuned from a Bose-Einstcin 
condensate to a state composed of localized atoms as the 
strength of a periodic optical potential was variedfii^ This 
is an example of a quantum phase transition; a phase 
transition generated by quantum fluctuations and corre- 
lations rather than by a competition between the energy 
of a system and the entropy of its thermal fluctuations 
Understanding this phenomenon has emerged as one of 
the most challenging and interesting tasks of condensed 
matter physics. Theoretically, it is generally accepted 
that it can be studied using the Bosc Hubbard model, 
where the transition is thought to be from a superfluid 
phase to a Mott-insulator (SF-MI), as examined in the 
seminal paper by Fisher et al.^ with an application to 
"^He absorbed in porous media in mind. The relevance of 
the Bose-Hubbard model to gases of alkali-metal atoms 
in optical lattices was realized in Ref. |5, and recent de- 
velopments have been reviewed in Refs. y and0- 

Interestingly, the Bose-Hubbard model is noninte- 
grable even in one dimension (as opposed to, say, its 
fermionic counterpart^). Gaining analytical insight into 
the SF-MI phase transition thus normally requires re- 
sorting to numerical and variational methods such as 
strong-coupling expansion;2ii£ coarse graining^ii mean- 
field theorieSfi^ field-theoretical approaches'^ or other 
perturbative methods'^ for a better understanding of this 
phenomenon. Within the variational approach, the phase 
transition is taken to be the point at which the variational 
ansatz has lower energy than a delocalized Bogoliubov 
state (where a fixed particle number at each lattice site 
is constrained). 

In a recent paper, Aizenman et al;^ considered an al- 
ternative model for the study of the SF-MI phase tran- 
sition. They studied the half- filled Bose-Hubbard model 



in the limit of infinite on-site repulsion (i.e., the case of 
hardcore bosons), in the presence of an alternating on- 
site chemical potential (a superlattice with period two). 
They showed that this model exhibits all the salient prop- 
erties apparent in the Bose-Hubbard model, while also 
being more 'treatable' analytically. Specifically, they 
were able to rigorously prove the existence of superfluid 
and Mott-insulating phases in three dimensions. In addi- 
tion, it is also known that this very same model is exactly 
solvable in one-dimension through a mapping to nonin- 
teracting fermions. In this case, the half-flUed system is 
insulating for any nonzero alternating potentialji^ The 
off-diagonal one-particle correlations and the momentum 
distribution function of this model, as well as its nonequi- 
librium dynamics, were computed by exact meansii in 
Ref. [H. 

Motivated by the aforementioned results, here we 
study the SF-MI phase transition of hardcore bosons in 
the presence of an alternating potential in two and three 
dimensions. We focus on the case where the system is 
at half-filling, in which case the transition between the 
superfluid state and the insulating state occurs at fixed 
density. Our first goal is to accurately determine the crit- 
ical values of the alternating potential strength at which 
the phase transition takes place. As a next step, we an- 
alyze the results of different approximate solutions, such 
as mean-field theory with and without the addition of 
spin-wave corrections, as these allow for an analytical 
treatment of the problem. 

Our approach is to first perform high-precision nu- 
merical simulations using the stochastic series expansion 
(SSE) algorith m^^i^° in order to find the critical points 
of the superfluid to Mott-insulator phase transition in 
the various dimensions. The quantities we calculate are 
the free energy $7, the density of bosons in the zero- 
momentum mode por^ the superfluid density ps and the 
compressibility k = dp/d^. The latter three quantities 
signify the transition from a superfluid to an insulator by 
dropping to zero at this point (while having nonzero val- 
ues in the superfluid regime) . We then employ mean-field 
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and spin- wave analyses, which allow for some analytical 
insight into the behavior of our observables of interest 
and the location of the critical point. Our use of these 
approximation methods is partly motivated by results 
previously reported by Bernardet et al.^ who studied 
the homogeneous version of the model in two dimensions. 
There, the mean-field approximation alone was shown to 
provide a fairly good qualitative description of the model, 
and remarkably enough, when spin- wave corrections were 
added, quantities such as the superfluid density and the 
condensate fraction were found to be virtually indistin- 
guishable from their exact-numerical counterparts. 

The paper is organized as follows. In Sec. [TI]we briefly 
review the model at hand. In Sec. IIIH we present the 
exact numerical solutions obtained using the stochastic 
series expansion (SSE) algorithm. Wc compute the vari- 
ous physical quantities at zero temperature, and find the 
critical values of the SF-MI phase transition. In Sec. 
IIVI we proceed to study several approximation schemes, 
namely mean-field approaches and spin-wave corrections, 
comparing the critical values obtained using these meth- 
ods, with the SSE results. In Sec. |V]we conclude with a 
few comments. 



II. THE MODEL 

The Hamiltonian for hardcore bosons in a period-two 
superlatticc in d-dimensions, with N = L"^ sites and pe- 
riodic boundary conditions, can be written as: 

(1) 

Here, (ij) denotes nearest neighbors, (al) destroys 
(creates) a hardcore boson on site fii — a\di is the local 
density operator, /i is the global chemical potential, and 
A(— l)'^^*) is an alternating local potential with a{i) = 
on the even sublattice and 1 on the odd sublattice. The 
hopping parameter t sets the energy scale. 

The hardcore boson creation and annihilation opera- 
tors satisfy the constraints 

af = a}^0, {a„al} = l, (2) 

which prohibit double or higher occupancy of lattice sites, 
as dictated by the — > oo limit of the Bosc-Hubbard 
model. For any two different sites i ^ j, the creation and 
annihilation operators obey the usual bosonic relations 

[a^, flj] = [a|, a]] = [a,, a]] = . (3) 

The expected phase diagram of the model, in dimen- 
sions higher than one, is sketched in Fig. [TJ Our model 
has two (trivial) insulating regimes corresponding to a 
completely filled lattice (with particle density p = 1), 
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FIG. 1: Qualitative description of the expected phase diagram 
of the model at hand, Eq. ([l]). The diagram contains three 
insulating regions corresponding to zero ('empty'), half ('MI') 
and full ('fully-filled') filling, and a superfluid (SF) phase. 



obtained for large and positive chemical potential values, 
and a second insulating regime which corresponds to an 
empty lattice, formed in the case where the chemical po- 
tential is large and negative. These two regimes are also 
present in the absence of the alternating potential. The 
alternating one-body potential creates another insidating 
phase, one for which the density of particles \s p = 1/2. 
In this case, the alternating potential, will in some cases 
(depending on its strength) create a gap in the energy 
spectrum, generating a superfluid to Mott-insulator tran- 
sition. As the latter regime is the one which is of interest 
to us, we shall henceforth set the global chemical poten- 
tial to /i = 0. In this case, the model has particle-hole 
symmetry which in turn fixes the density at p = 1/2 as 
desired. 

Before moving on, a remark is in order. The p = 1/2 
insulating phase of the model at hand is a consequence 
of a counterbalance between strong on-site interactions 
(which in our model arc in fact infinite) and an alter- 
nating potential. The resulting local density will thus 
be different on the odd sublattice than on the even sub- 
lattice. While this state is sometimes referred to as a 
charge density wave,— in what follows, we shall address 
this phase as a Mott-insulator, in the spirit of Rcf. ITsl 



III. NUMERICAL RESULTS 

We obtain numerically-exact results for the model at 
hand by performing numerical simulations based on the 
stochastic series expansion (SSE) algorithmJ^i^ As our 
main objective is to find the critical points of the SF-MI 
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phase transition in the various dimensions, simulations 
are performed for a range of values of the ratio A/t and 
for various system sizes. Since we arc interested in the 
zero-temperature properties of the system, simulations 
are performed with high inverse-temperature (3 = l/T 
(in our units, fcs = 1), where in most cases we will find 
it sufficient to have /3 > 2L in order to obtain virtually 
zero-temperature results. (The effects of increasing (3 
beyond this value are indiscernible.) 

Finite size effects are corrected by repeating the simu- 
lations with different system sizes. The thermodynamic- 
limit value of the phase transition is then extrapolated by 
performing finite size scaling of the results in the vicinity 
of the phase transition: Around the critical point, most 
physical quantities (which we denote here by X) scale 
according to the general rule: 

XL^/" = F{\A~ A^\L^/''), (4) 

where is a universal scaling function, A — Ac is the 
shifted control parameter (A being the control parame- 
ter, and Ac - the critical value), v is the correlation length 
critical exponent and ^ is the critical exponent belonging 
to the observable X. The values of these exponents are 
determined by the universality class the transition be- 
longs to. In our case (and in the Bose-Hubbard model 
for integer filling as well), it is the (d+l) dimensional XY 
universality class4i^ We note that the above universality 
class characterizes only the fixed-density transition (the 
dashed line in Fig.[T]). The transition driven by changing 
the density belongs to the mean-field universality class 
and is characterized by different critical exponents. 

Equation Q above will help us find the critical point, 
as it tells us that (a) the quantity XL^I'^ should be inde- 
pendent of the system-size at the phase transition, and 
(b) when plotting XL^I'' against \A - Ac\L^I" the re- 
sulting curve should be independent of the system-size 
as well. 

The quantity we shall be using to that end is the su- 
pcrfluid density, which has the critical exponent ^ — 
v(d + 2 — 2) (see Ref. for details) where d is the di- 
mension, and z is the dynamical critical exponent, which 
in our case is z = li^ The correlation length exponent 
V is dimension-dependent and takes the values 1, 0.672 
and 0.5 in one, two and three dimensions, respectively. 

A. One dimension 

In one dimension, our model has an analytic solutioUiiS, 
This is due to the Jordan- Wigner transformation which 
enables the mapping of the hardcore bosons Hamiltonian 
to that of noninteracting spinless fermionsii^ The latter 
Hamiltonian may be diagonalized to produce exact ana- 
lytical results. In this case, the SF-MI phase transition 
occurs at Acji^dC) = 0, i.e., the system is superfluid only 
when the alternating potential is absent, in which case 
it exhibits off-diagonal quasi-long-rangc order (a power- 
law decay of the onc-particle correlations). In that sense. 



one may say that the system exhibits quasi-condensation 
when A = OMOS 

Simulations in one dimension were thus performed only 
as a check on our computational method. No discrep- 
ancies between the analytical solution and the numeri- 
cal one were found: In Fig. [2l the superfluid density is 
plotted against Aji^dt) for different system sizes (here, 
(i = 500). In the figure, all curves intersect at the critical 
point Acji^dt) = 0, indicating the location of the phase- 
transition in the thermodynamic limit, in agreement with 
the analytic results. The inset shows the scaled superfluid 
density as a function of the scaled control parameter, in 
which case all curves should be, and in fact are, on top of 
each other. The numerical value for the superfluid den- 
sity at the transition coincides with the expected value 
of 7r~^ given by the analytic solution^i^ 
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FIG. 2: (Color online) Scaled superfluid density as a function 
of A/(2dt) for the various system sizes in the one-dimensional 
case. The intersection at A/{2dt) — indicates the location 
of the SF-MI phase transition. In the inset, the control pa- 
rameter (the horizontal axis) is scaled as well, leading to the 
collapse of all data points into a single curve. 

As superlattices such as the one we study here have al- 
ready been realized in experiments with ultracold bosons 
in optical lattices,— i^i2^>^ and the observable usually 
measured in those kind of experiments is the momentum 
distribution function n{k), wc plot it in Fig.[3]for two dif- 
ferent values oi A/t. Due to the quasi- long-range decay of 
one-particle correlations in the superfluid phase, the mo- 
mentum distribution function has a peak at fc = [Fig. 
[31[a)]. On the other hand, in the insulating phase, the 
one-particle correlations decay exponentially, yielding a 
broad momentum distribution [Fig. [3]Jb)]. This leads to 
the following observation: As one increases the size of the 
lattice (while keeping the density fixed), one finds that 
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in the superfluid phase n{k) increases for small values of 
k [Fig. Elja)], while for the insulating phase [Fig. ^h)] 
this docs not happen. Exact results for n{k) (using the 
approach described in Ref . ITtI) . are also presented in Fig. 
[3l As expected, the SSE results are right on top of the 
exact ones. 




FIG. 3: (Color online) Momentum distribution function n{k) 
in the superfluid regime (top) and in the insulating regime 
(bottom) for the one-dimensional system with 100 sites. In 
one dimension, the system is superfluid only at A = (top 
panel). This is shown by the sharp peak in the k = mode 
of the momentum distribution function which diverges in the 
thermodynamic limit. In this case, the system exhibits quasi- 
long-range order. In both panels, the SSE results (empty 
circles) are on top of the analytical ones (full circles), serving 
as an indication to the accuracy of our computational method. 



The value for the critical point we obtained here agrees 
with the value recently obtained by Priyadarshee et ali^ 
The momentum distribution function in the superfluid 
and insulating regimes are shown in Figs.[5lja) and[5l^b), 
respectively. In two dimensions, the superfluid regime 
exhibits true off-diagonal long-range order, which means 
that the n{k = 0) peak is sharper that in one dimen- 
sion, which exhibits only quasi-long-range order. This 
can be seen in Fig. El^a). The Mott-insulating regime 
is once again characterized by an exponential decay of 
one-particle correlations. The corresponding momentum 
distribution function has a broad peak around n{k = 0) 
as shown in Fig. [^b). 
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FIG. 4: (Color online) Scaled superfluid density as a function 
of A/{2dt) for the various system sizes in the two-dimensional 
case. The intersection at Ac/{2dt) « 0.495 indicates the oc- 
currence of the phase transition at that point. In the inset, 
the control parameter (the horizontal axis) is scaled as well, 
leading to the collapse of all data points into a single curve. 



B. Two dimensions 

In dimensions higher than one, no analytic solution 
to the model exists, so accurate results are obtainable 
only numerically. Here, we have applied the SSE al- 
gorithm to systems of sizes ranging from 10 x 10 to 
64 X 64, with inverse-temperature f3 — 96. In Fig. IH 
the scaled superfluid density is plotted against A/{2dt) 
for the different system sizes (the errors are on the order 
of magnitude of the symbol sizes). All curves intersect 
at Ac/{2dt) = 0.495(±0.004), signifying the phase tran- 
sition. The inset shows the scaled superfluid density as a 
function of the scaled control parameter. As in the one- 
dimensional case, all data points fall into a single curve. 



C. Three dimensions 

In three dimensions, we have performed simulations 
for system sizes ranging from 6 x 6 x 6 to 20 x 20 x 20 
and an inverse temperature of /3 = 40. Here, the SF-MI 
phase transition is found at Ac/{2dt) = 0.693(±0.005), 
as indicated by the scaled superfluid density plotted as 
a function of A/{2dt) in Fig. [SI for the different system 
sizes. The inset in Fig. [S] depicts the scaled superfluid 
density as a function of the scaled control parameter, 
exhibiting the collapse of all data points into a single 
curve, as in one and two dimensions. The momentum 
distribution function in three dimensions is qualitatively 
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FIG. 5: (Color online) Momentum distribution function n{k) 
in the superfluid regime A/{2dt) = 0.1 (top) and in the in- 
sulating regime A/{2dt) = 0.7 (bottom) for a 64 x 64 system 
and f3 = 96. 
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FIG. 6: (Color online) Scaled superfluid density as a func- 
tion of A/{2dt) for the various system sizes in the three- 
dimensional case. The intersection at A^/ {2dt) ~ 0.693 in- 
dicates the location of the SF-MI phase transition. In the 
inset, the control parameter (the horizontal axis) is scaled as 
well, leading to the collapse of all data points into a single 
curve. 



alternating magnetic field applied along the z direction: 



similar to its two-dimensional counterpart, both in the 
superfluid phase and in the insulating phase, and thus 
will not be presented here. 



IV. APPROXIMATION SCHEMES 



H = - tJ2{S+Sr +S+S- 



5][/.+^(-ir« 



(6) 



Having obtained the critical values via quantum Monte 
Carlo techniques, we now turn to look for approxima- 
tion schemes that would provide analytical descriptions 
of the phase transition. We start this investigation with 
the Gutz wilier mean-field approach. Before doing so 
however, we recall that the model at hand can also be 
viewed as the XY model of a spin-1/2 system<^ We shall 
make use of this correspondence, utilizing the exact map- 
ping between bosonic operators and SU{2) generators, 
namely, 

4 St, (5) 

With this mapping, the hardcore bosons Hamiltonian, 
Eq. ID), becomes that of the XY antiferromagnet with an 



A. Mean-field approach 

We start our mean-field calculation with the following 
product state as an initial ansatz: 



where {9j,(pj) specify the orientation of the j-th spin. 
Obviously, we expect every other site to be described 
by the same wave function, due to the symmetry of the 
problem. This is schematically shown in Fig. [T] As we 
are using the grand-canonical scheme, the orientations of 
the spins will be determined by minimizing the grand- 
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canonical potential (per site) 



MF 



(0|iJ|0) 



MF 



sm til sint/w cos 



-±j2[^+A{-ir^'^^] (1+COS0O 



0.) 

(8) 



with respect to these angles. For the azimuthal an- 
gles, this simply implies a constant (yet arbitrary) value 
ifj = while for the polar angles, extremization yields 



COS^i 



cos O2 




(9a) 
(9b) 



where ^1^2 = A)/{2dt). These values correspond to 
a minimal configuration only in the region |/ii/X2| < 1. 
Outside this region, the system is saturated, and the so- 
lution is one where all spins are aligned - pointing ei- 
ther all up or all down. In bosonic language, these latter 
configurations correspond to the completely full/empty 
insulating configurations. 

At this point we can calculate the following quantities. 
First, the density of particles is: 



Pmf 



N 

= i + i (cos 6*1 + cos6'2) 



(10) 



Next, the free energy becomes 

r^MF 



MF(0|i?|0)MF = sin 6*1 sin 6*2 - ^ 



2 

- i (^ + A)cos6'i - i (/i - ^)cos6'2 , (11) 
and the density of bosons in the zero-momentum mode 
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FIG. 7: A schematic description of the product state in the 
mean-field approach in two dimensions. Every other site is 
described by the same wave function. 



po is calculated as: 
1 



PO.MF 



N 



MF 



(0|aj^^gafc=o|0) 



MF 



(12) 



4iV 



-y 



sm f i sm ( 



1 

16 



(sin 9i 



sm O2 ) 



The superfluid density ps requires a special treatment 
of the boundary conditions. As is well knowu)^ one can 
relate the superfluid density to the "spin stiffness" . To 
accomplish this, one needs to compare Vt (the free en- 
ergy) of the system under periodic conditions with the 
free energy under a "twist" in the boundary conditions 
along one of the linear directions (say, the x direction). In 
the periodic case, which we treated above, the azimuthal 
angles Lpj were all identical. To implement a twist, we 
take this angle to be site-dependent and with a constant 
gradient such that the total twist across the system in 
the X direction is tt, namely 5lp = ^Pj+x ~ = tt/L. 
Within the mean-field treatment, one can show that 
addition of this gradient is equivalent to substituting 
t ^ t/d[{d — 1) + cos 6cp]. Now, the square of the gradi- 
ent twist is related to the superfluid density viaj^ii^ 



'twisted — = tpsStp'^ , 

which in turn yields the simple expression 

1 dfl 



(13) 



(14) 



Setting /i = 0, this expression for the superfluid density 
coincides with that of Pq^mf- 



Ps,MF = PO,MF — 







A Y 



•.Adt 



< 1 



A 

^>1 
2dt — ^ 



(15) 



giving the critical value for the phase transition 
Ac/{2dt) = 1. Figures [5] and [5] show: (a) the free en- 
ergy, (b) the superfluid density, (c) the density of bosons 
in the zero-momentum mode, and (d) the compressibil- 
ity of the system as a function of A/{2dt) in two and 
three dimensions. The dashed and solid curves represent 
the mean-field and SSE results, respectively. As one can 
immediately see, the critical values obtained within the 
mean-field approximation do not agree with the exact- 
numerical results. In two dimensions the error is ~ 100% 
and in three dimensions, it is w 50%. The very large 
errors here merely reflect the fact that the mean-field ap- 
proach used here is not fit to describe the model at hand, 
especially in the vicinity of the SF-MI phase transition. 



B. Adding spin- wave corrections 

As pointed out earlier, the addition of spin- wave cor- 
rections yields virtually exact results in the homogeneous 
case in two dimensionsi^ For the reader's convenience, 
we review the mean-field calculations of the homogeneous 
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FIG. 8: (Color online) Thermodynamic quantities in two di- 
mensions, (a) Free energy [t — 1], (b) superfluid density, (c) 
the density of bosons in the zero-momentum mode, and (d) 
compressibility as a function of A/{2dt). The solid lines in- 
dicate the SSE results (64 x 64 sites, /3 = 96), whereas the 
dashed, dotted and dash-dotted lines indicate the mean-field, 
mean-field plus spin-waves and improved mean-field results, 
respectively. 



FIG. 9: (Color online) Thermodynamic quantities in three di- 
mensions, (a) Free energy [t — 1], (b) superfluid density, (c) 
the density of bosons in the zero-momentum mode, and (d) 
compressibility as a function of A/{2dt). The solid lines indi- 
cate the SSE results (16 x 16 x 16 sites, /3 = 40), whereas the 
dashed, dotted and dash-dotted lines indicate the mean-field, 
mean-field plus spin-waves and improved mean-field results, 
respectively. 



(A = 0) case and its spin-wave corrections in Appendix 
El (thereby also correcting some misprints that appeared 
in the original manuscript examining this case, Ref. [2lh . 
Let us see how the mean-field results are modified by the 
addition of spin-wave corrections in our case. To include 
these, we proceed in the usual wayj^'^i^'^^ We first in- 
troduce a set of local rotations that align the z direction 
of each of the spins with its mean field orientation. This 
is accomplished by switching to new spin operators de- 
fined by 



■5f 

S'y I ^R{9,,^,) 
ST 




(16) 



where R{Oj, ^Pj) is the 3x3 rotation matrix 



' cos dj COS Lpj 

COS Q j sin ipj cos 
— sin Qj 



— sin Lpj sin Qj cos ipj 
cos sin Qj sin (pj 
cos Qj 








The corresponding new annihilation and creation oper- 
ators hj ^ S'~ and ^ describe low-energy fluc- 
tuations about the mean-field ground state - these are 
spin waves. They too obey hardcore bosons commuta- 
tion relations. Substituting these expressions into our 
Hamiltonian, and ignoring cubic and quartic terms in 



these bosonic operators (thus assuming a dilute gas of 
spin waves), the new Hamiltonian reads 

^sw = i?MF + ^^6tS, + C^(-l)^«6^6, (18) 

i i 



(ij) 

where the coefficients are 



A = t(l + cos6'icos6'2) , (19a) 

B = f/2(l-cos6'iCos6'2) , (19b) 

C ~ dt {fj,i cos 9i — ^2 cos 62) , (19c) 
D = (it (2 sin6'i sin + Ml cos 01 4- cos ^2) ■ (19d) 

This quadratic Hamiltonian can be diagonalizcd by first 
going to Fourier space, using bj ~ N~^/^ J^A: 
This in turn yields the Hamiltonian: 

Hsw^Hmf + Y.{D~Ajk)blbk + CY,blh+L/2 

k k 

+ BY,ik{bibi_,+bkbL^k), (20) 



where, 7^ = J2i=i cos (^x^) , and ki . . . kd are the com- 
ponents of the momentum vector in each of the di- 
rections. We note that the Fourier-space operators bk 
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and ol no longer obey the hardcore bosons commuta- 
tion relations. These field operators are only excita- 
tions about the ground state, and are treated as soft-core 
bosons i^ii^i^i^i^ At this point, our Hamiltonian may 
be diagonalized in a straightforward manner (we review 
the diagonalization process in Appendix |B| . Once diag- 
onalized, the Hamiltonian takes the form 



^sw = ^^MF + ^ A 



En 



(21) 



where the A^ 's are energy levels and Eq is the correction 
to the ground-state energy of the system, given by: 



E. = -J2^-2D 



-f- \ (A2 ~ 4B2)^2 ^ £,2 _^ _ 2J{DCY + [{ADY - {2BCf]^l 



(22) 



The operators fy^ and fjk in Eq. (|2ip arc mod- 
ified spin-wave creation and annihilation operators, 
respectively, and are each a linear combination of 
bk,bL-k,bk+L/2,bL/2-k and their adjoints. The coeffi- 
cients of these linear combinations are fixed during the 
diagonalization process, and using them, all physical ob- 
servables can be calculated in a straightforward manner 
(we elaborate on this matter in Appendix IB|) . 

The results of the spin-wave analysis are indicated 
by the dotted lines in Figs. [5] (two dimensions) and [H] 
(three dimensions). They show: (a) the free energy, (b) 
the superfluid density, (c) the density of bosons in the 
zero-momentum mode, and (d) the compressibility, after 
the addition of spin-wave corrections, as a function of 
A/{2dt). 

As one can see in those figures, in the superfluid 
phase, the spin-wave corrected values for the free en- 
ergy are almost on top of the exact-numerical ones; and 
more so in the three-dimensional case than in the two- 
dimensional one. As for the other measured observables, 
the spin-wave corrections are clearly an improvement 
over the mean-field results, especially for small values of 
A/t where the spin- wave corrections yield virtually ex- 
act results. Unfortunately however, as one approaches 
the phase transition itself, the spin-wave corrections lose 
their accuracy, eventually leaving the phase-transition at 
its mean-field value, namely at Ac/{2dt) = 1. 

Another issue worth noting here is the behavior of the 
spin- wave corrected superfluid density [Figs. [SJ^b) and 
Wijo)] in the vicinity of the predicted phase transition, 
A/{2dt) ~ 1. On the superfluid side of the transition 
the superfluid density becomes negative, indicating the 
breakdown of the spin-wave approximation for that quan- 
tity. The transition point is still signaled by a discontinu- 
ity in ps- However, the overall behavior of the superfluid 
density around the transition point is clearly an artifact 
of the spin-wave approximation and should not be con- 



sidered further. 



C. Improved mean-field approach 

Having seen that spin-wave corrections, albeit accurate 
in the weak-potential regime, do not modify the critical 
point predicted by the mean-field solution, we have de- 
vised an improved mean-field approach. As wc show now, 
this method provides a significant improvement over the 
mean-field results (and the spin-wave corrections) dis- 
cussed previously, particularly in the context of the loca- 
tion of the critical point. 

We start with a variational ansatz which, as before, is 
a product state. However, this time we do not choose a 
product of single-site wave-functions. The new ansatz is 
a product of wave- functions each describing the state of 
a 'block' of 2*^ sites, such that with this block as the basic 
cell, the model turns homogeneous. In two dimensions a 
block consists of 2 x 2 cells (as shown in Fig. [TO)) each of 
which is described by the general wave function 

|0)iMF= n ( (^^Jkl\ijkl)] , (23) 

blocks \ij,fc,;6{i,T} / 

where the generalization to three dimensions, in which 
case the basic block is a 2 x 2 x 2 cell, is straightfor- 
ward (note that the coefficients for each of the blocks 
are the same). As before, we minimize the free energy 
^^IMF = iMF(0|-ff |0)iMF with respect to the coefficients 
Cijki of the wave function (this time we do so numer- 
ically). Obtaining the various observables in terms of 
the wave function given in Eq. ((23l) is straightforward, 
and was performed in much the same way as the usual 
mean-field approach discussed in Sec. IIV Al The results 
of this approximation arc given by the dash-dotted lines 
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in Figs. [8] and [9l They depict: (a) the free energy, (b) 
the superfluid density, (c) the density of bosons in the 
zero-momentum mode, and (d) the compressibility, as a 
function of A/{2dt). 

As the figures indicate, in most instances, the results of 
this method are more accurate than those of the previous 
approximation schemes, in particular, for the location of 
the phase transition. The critical values given by this 
approximation arc Ac/{2dt) ~ 0.815 in two dimensions 
{^i 60% error) and Ac/{2dt) = 0.875 in three dimensions 
(« 24% error). Also, we note that while the spin- wave 
corrected values for the various thermodynamic quan- 
tities are a better approximation in the weak potential 
[small A/{2dt)] regime, as one moves away from this re- 
gion, the improved mean-field technique proves to be a 
better estimator for all quantities but the free energy. It 
is clear still that the improved-mcan-field method pre- 
sented here is far from being very accurate. 







1 





FIG. 10: In the 'improved mean-field' larger unit-cell 

is defined. In the two-dimensional case at hand, the new cell 
consists of 2 X 2 sites. With this new definition the model 
turns homogeneous and a product of identical wave functions 
is then guessed as a solution. 



V. CONCLUSIONS 

We have studied the superfluid to Mott-insulator phase 
transition of hardcore bosons in a period-two superlatticc 
in two and three dimensions. We focused on the case 
where the system is at half filling, for which the quantum 
phase transition belongs to the {d + 1) dimensional XY 
universality class. 

Using quantum Monte Carlo simulations and finite size 
scaling, we have determined the critical value of the alter- 
nating potential parameter A at which the SF-MI phase 
transition occurs. In two dimensions, our results agree 
with previous calculations 

We have also compared our numerical results against 
several approximation schemes, some of which have been 
used successfully in the two-dimensional homogeneous 
version of the model. We have seen that employing a 
mean-field approach using the usual Gutzwiller ansatz 
works very poorly (w 100% error in two dimensions and 



w 50% error in three dimensions). This is a clear indi- 
cation of the fact that this mean-field approach is not 
suitable for describing this model, especially in the vicin- 
ity of the phase transition, as it breaks down in the strong 
coupling regime. 

The spin-wave corrections to the mean-field solution 
turned out to be very useful, especially in the superfluid 
phase, where the spin-wave corrected estimation of the 
free-energy is very close to the exact values, and also 
reproduced the exact results for all observables for small 
values of A/t. However, as one moves away from the weak 
potential regime, the spin- wave corrections become more 
and more inaccurate, and their predictions of the critical 
points eventually coincide with those of the mean-field 
approach, therefore indicating their unusefulness in that 
region. 

The improved mean-field approximation scheme we 
have devised here, which was based on the underlying 
homogeneity of the problem, has proved to be an im- 
provement over the previous methods, albeit still far 
from being accurate. This approach provides an analyt- 
ical description of the superfiuid-to-Mott-insulator tran- 
sition and gives an estimate of the critical value for the 
transition with (around) one half the error of the usual 
Gutzwiller ansatz, i.e., it is an improvement in terms of 
the location of the critical point. 
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APPENDIX A: MEAN-FIELD RESULTS AND 
SPIN- WAVE CORRECTIONS IN THE 
HOMOGENEOUS CASE 

In what follows, we briefly review the results of the 
mean- field approximation of Sec.[lVX]and its spin- wave 
corrections in the homogeneous {A = 0) case for arbitrary 
values of ^. 

Starting with the ansatz given in Eq. ([7]), minimiza- 
tion of the free energy Eq. ([8|) with respect to the spin 
orientation angles yields 



cos 6-i 



2dt 



(Al) 



where the azimuthal angle takes on, once again, a con- 
stant yet arbitrary value (/9j = The density of particles 
becomes 



N 

2^^1dt) ' 



(A2) 
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and the free energy is 



MF 



MF 



{0\H\0) 



MF 



dt , ^ I , 

— sm fc*— -^(l + cos( 



— V 

2dtJ 



(A3) 



The density of bosons in the zero-momentum mode turns 
out to be 



PO,MF — -^MF(0|aj^^gafc=o|0)MF 

1 

4 



(A4) 



4iV2 



sm f , sm ( 



\2dtJ 



Using Eq. (|14p . it can be easily shown that the expression 
for the superfluid density Ps.MF in the homogeneous case 
coincides with the expression obtained for po,MF above 
(as is the case with the alternating potential). 

The addition of spin- wave corrections to the mean-field 
results is carried out in exactly the same manner as with 
the staggered potential. The Hamiltonian in this case 
has the same form as the one given in Eq. but with 
coefficients 



A = t 


1 + 


f— ) 

\2dtJ 


(A5a) 


B^i 

2 


1 - 


-(—Y 

\2dtJ 


(A5b) 


C = 0, 






(A5c) 


D = 2dt. 




(A5d) 



The spin-wave field operators which diagonalize the 
Hamiltonian are given by the simple relation)^ 



Sfe = cosh (j)k 7% - sinh cfikVL^k' 
with (pk obeying 



sinh 



cosh 



2 \ ^{D-A-fkY-[2B^k? 

b = - ( D-Alk 
' 2 [^iD-Any-{2B^,y 



(A6) 

(A7a) 
(A7b) 



and so, the various spin-wave corrected quantities may 
be written explicitly: the corrected density of particles is 



Psw = Pmf 



liLysinh^ 
N2dt ^ 



(A8) 



and the free energy becomes 

^sw = ^MF (A9) 



Using Eq. (fT4|) . the superfluid density immediately fol- 
lows. Finally, the density of bosons in the zero- 
momentum mode turns out to be: 



1 

PO,SW — PO.MF ^ 



1 - 



(—Y 

\2dt) 



sinh^ (j)k ■ 



(AlO) 



APPENDIX B: DIAGONALIZATION OF 
QUADRATIC BOSONIC HAMILTONIANS 

Following the prescription given in Ref. [13 for the di- 
agonalization of fermionic quadratic Hamiltonians, we 
provide here the analogous prescription for the diagonal- 
ization of bosonic quadratic Hamiltonians of the general 
form 



k.rn 



(Bl) 



where bk and b\ are bosonic annihilation and creation 
operators, respectively, and Ak„i and Bk„i are real- valued 
and symmetric. For the spin- wave Hamiltonian of Eq. 
([20|) . the coefficients are 

Akrn = {D - A^k)5km + C4_,„+i/2 , (B2a) 
Bkra = 2B-ik5k.L~ra ■ (B2b) 

The diagonalization process starts by defining the follow- 
ing linear transformation 



(B3) 



where gkm and hkm are real-valued and we require f/k 
and fij. be bosonic operators. This is enforced by the 
constraint 

Skm = [Vk,1lL] = E idkig-ml - hklhml) ■ (B4) 



The coefficients gkm and hj^m are determined in such a 
way that the transformed Hamiltonian takes the diagonal 
form 



H 



E 



(B5) 



once the fjkS are substituted for the 5fe's, and 
Eq = "Efcm^Ti'^mfc- ^lic ncw Hamiltonian is al- 
ready in diagonal form, the new field operators obey the 
eigenvalue equation 

[fjk,H]=Kkfik. (B6) 

Plugging in the transformations given in Eqs. (jBSp . we 
obtain the relations 

Afefffcrn = ^ {gkiAirn - hklBml) , (B7a) 
/ 

^khkm — ^ {gklBml — hklAim) ■ (B7b) 
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These relations may be further simpHfied by defining the 
new coefficients 

^fcm = gkm + hkm , (B8a) 
Ipkni = gkm - hkm , (B8b) 

for which, the constraint (|B4p translates to 

^ X! ('^kli'ml + i'kl<l)ml) = Skm ■ (B9) 

With the above definitions, Eqs. (jB7p may be cast in 
vector notation: 

0,(A-i?)=Afet/>,, (BlOa) 
V,(A + B)=Afe</.,. (BlOb) 

These can be solved by simply plugging each of these 
equations into the other, resulting in the eigenvalue equa- 
tions 

^Pk{A + B){A-B) =AliPk, (Blla) 
ct>k{A-B){A + B) =Al<p,. (Bllb) 



These equations are to be solved by standard techniques. 
Once the A^'s, (p^s and ip^s are found, all physical ob- 
servables can be readily calculated: First, the observ- 
able of interest should be expressed in terms of normal- 
ordered %'s. This may be accomplished by using the 
inverse of the transformation given in Eq. (|B3p : 

(B12) 

As a next step, one should use the fact that as excitations, 
the 77fc's obey 77a;|0)mf = 0. This leads to 

I 

(B13) 

As an example, consider the spin-wave corrected density 
of particles in our model. It is calculated as 



PSW = ^ MF(0|alat|0)MF = PMF - ^ MF{0\blk COs6>a|0)MF (E 



PMF - ^ (cos 6*1 -f COs6'2)^MF(0|6[.6fe|0)MF " ^ ( COS 6*1 - C0S6'2) ^MF(0|6[.6fc+i/2|0)MF 

fc k 

- ^ I (cos 01 + COS 02) J2 i^km ' i'km)^ + (^OS ^1 " COS02) ('^fcm " '^fc™) {'t'ik+L/2).m " "(-'{k+Lm, 

\ m,k mk 



All other observables may be calculated in the same manner. 
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